/*-------------------------------------------------------------------------*/
/**
  @file		match.c
  @author	Renaud Wanchoor
  @date	July 2004
  @version	$Revision: 2.13 $
  @brief	matching schemes

  Implementation of various matching schemes
  all prototypes should be the same
  \param coactrl coarsening control structure
  \return coarsened graph

*/
/*--------------------------------------------------------------------------*/
#include "pagrid.h"


static int __Match_LAM_opt_w3 (CoaCtrl* coactrl, int* red_n);
static int __Match_LAM (CoaCtrl* coactrl,
			int a, int b, int edgewab, 
			int* edge_p_check, 
			int *red_n);

/**
 * \brief LAM matching scheme
 *
 * LAM Implementation, stolen from PARTY
 * \author Robert Preis, Universit\"at Paderborn, Germany
 * preis@hni.uni-paderborn.de
 *
 * modified by Renaud WANSCHOOR 2004/1/29 to use PAGRID data structure and
 * random numbers
 * \param coactrl coarsening control structure
 * \return coarsened graph 
 */
GraphType* LAM(CoaCtrl* coactrl)
{
  GraphType *graph = coactrl->g;

  int
    cnvtxs         = 0,
    randguide      = coactrl->randguide,
    n              = graph->nvtxs,
    red_n          = graph->nvtxs,
    *edge_p        = graph->adjidx,
    *edge          = graph->adjncy,
    *edge_w        = graph->adjwgt,
    *matching      = graph->match,
    *cmap          = graph->map,
    *edge_p_check  = U_malloc(n+1,int,"Match_LAM error code 1\n");
  
  
  int 
    ii, i, j, 
    *perm          = coactrl->guide; 

  for (i=0; i<n; i++)  
    matching[i] = i;

  memcpy(edge_p_check,edge_p,(n+1)*sizeof(int));

  if(randguide == 1)
    randnum(n, perm);
  else
    for(i=0;i<n;i++)
      perm[i] = i;
  
  for (ii=0; ii<n; ii++){

    i = perm[ii];

    for (j=edge_p[i]; matching[i]==i && j<edge_p[i+1]; j++){
      if (matching[edge[j]] == edge[j]){
	__Match_LAM(coactrl,
		    i,edge[j],edge_w[j],
		    edge_p_check,&red_n);
      }
    }
  }

  __Match_LAM_opt_w3(coactrl,&red_n);

  /* build cmap */
  for (i=0;i<n;i++)
    cmap[i] = -1;

  for (ii=0; ii<n; ii++){
    i = perm[ii];
    if(cmap[i] != -1)
      continue;
    cmap[i] = cmap[matching[i]] = cnvtxs++;    
  }

  graph = Contract(graph, perm, cnvtxs);

  U_free("Match_LAM error code 2\n",&edge_p_check,LEND);

  return graph;
}
/** \brief lam recursive function */
int __Match_LAM (CoaCtrl* coactrl,
		 int a, int b, int edgewab, 
		 int* edge_p_check, 
		 int *red_n){
  GraphType *graph = coactrl->g;

  int
    redl           = coactrl->npro,
    *edge_p   	   = graph->adjidx,
    *edge          = graph->adjncy,
    *edge_w        = graph->adjwgt,
    *matching      = graph->match;
  
  int	
    edge_p_check_a=edge_p_check[a],
    edge_p_check_b=edge_p_check[b],
    neigh, neigh_w, neigh_deg, 
    edge_p_a1=edge_p[a+1], 
    edge_p_b1=edge_p[b+1],
    deg_a=edge_p[a+1]-edge_p[a], 
    deg_b=edge_p[b+1]-edge_p[b];

  while (matching[a]==a && matching[b]==b &&
	 (edge_p_check[a]<edge_p_a1 || edge_p_check[b]<edge_p_b1) &&
	 (*red_n)>redl ){
    if (edge_p_check[a]<edge_p_a1){
      neigh     = edge[edge_p_check[a]];
      neigh_w   = edge_w[edge_p_check[a]];
      neigh_deg = edge_p[neigh+1]-edge_p[neigh];
      (edge_p_check[a])++;
#ifdef LAM_ORIG
      if (neigh_w>edgewab)
#else
      if ((deg_b>1&&neigh_deg>1&&(neigh_w>edgewab||(neigh_w==edgewab&&neigh_deg<deg_b))) ||
	  (deg_b==1&&neigh_deg==1&&neigh_w>edgewab) ||
	  (deg_b==1&&neigh_deg>1&&neigh_w>2*edgewab) ||
	  (deg_b>1&&neigh_deg==1&&2*neigh_w>=edgewab) )
#endif
        __Match_LAM(coactrl,
		    a,neigh,neigh_w,
		    edge_p_check,red_n);
    }
    
    if (matching[b]==b && edge_p_check[b]<edge_p_b1){
      neigh     = edge[edge_p_check[b]];
      neigh_w   = (edge_w?edge_w[edge_p_check[b]]:1);
      neigh_deg = edge_p[neigh+1]-edge_p[neigh];
      (edge_p_check[b])++;
#ifdef LAM_ORIG
      if (neigh_w>edgewab)
#else
      if ((deg_a>1&&neigh_deg>1&&(neigh_w>edgewab||(neigh_w==edgewab&&neigh_deg<deg_a))) ||
	  (deg_a==1&&neigh_deg==1&&neigh_w>edgewab) ||
	  (deg_a==1&&neigh_deg>1&&neigh_w>2*edgewab) ||
	  (deg_a>1&&neigh_deg==1&&2*neigh_w>=edgewab) )
#endif
	__Match_LAM(coactrl,
		    b,neigh,neigh_w,
		    edge_p_check,red_n);
    }
  }

  if (matching[a]==a && matching[b]==b && (*red_n)>redl){
    matching[a] = b;
    matching[b] = a;
    (*red_n)--;
  }else if (matching[a]==a && matching[b]!=b)
    edge_p_check[a] = edge_p_check_a;
  else if (matching[a]!=a && matching[b]==b)
    edge_p_check[b] = edge_p_check_b;

  return 0;
}
/** \brief lam optimization function (augmenting path of length 3) */
int __Match_LAM_opt_w3 (CoaCtrl* coactrl, int* red_n){
  int	
    i, j, k, 
    free_p=0, 
    matching_weight,
    gain_2, best_2=-1, 
    best_2_w, gain_3, 
    best_near=-1, 
    best_middle=-1,
    best_distant=-1, 
    match_edgew_ab, match_edgew_cd;

  GraphType *graph = coactrl->g;

  int 
    redl        = coactrl->npro,
    n           = graph->nvtxs,
    *edge_p     = graph->adjidx,
    *edge       = graph->adjncy,
    *edge_w     = graph->adjwgt,
    *matching   = graph->match,
    *stack_free = graph->map;
  

  for (i=0; i<n; i++) 
    if (matching[i] == i)
      stack_free[free_p++] = i;
  
  while (free_p && (*red_n)>redl){
    i = stack_free[--free_p];
    gain_2 = gain_3 = -1;
    for (j=edge_p[i]; matching[i]==i && j<edge_p[i+1]; j++){
      if (matching[edge[j]]==edge[j]){
	matching[i] = edge[j];
	matching[edge[j]] = i;
	j=edge_p[i+1];
	(*red_n)--;
      }else{
	for (k=edge_p[edge[j]]; matching[edge[j]]!=edge[k]; k++) ;
	matching_weight = edge_w?edge_w[k]:1;
	if ((edge_w?edge_w[j]:1)-matching_weight > gain_2){
	  gain_2 = (edge_w?edge_w[j]:1)-matching_weight;
          best_2 = edge[j];
          best_2_w = edge_w?edge_w[j]:1;
        }
	  
        for (k=edge_p[matching[edge[j]]]; k<edge_p[matching[edge[j]]+1]; k++)
          if (matching[edge[k]]==edge[k] && edge[k]!=i)
            if((edge_w?edge_w[j]:1)-matching_weight+(edge_w?edge_w[k]:1)>gain_3){
	      best_near = edge[j];
	      best_middle = matching[edge[j]];
	      best_distant = edge[k];
	      gain_3 = (edge_w?edge_w[j]:1)-matching_weight+(edge_w?edge_w[k]:1);
              match_edgew_ab=edge_w?edge_w[j]:1;
	      match_edgew_cd=edge_w?edge_w[k]:1; 
	    }
      }     
    }

    if (matching[i]==i){
      if (gain_3 >= 0){
	matching[i] = best_near;
	matching[best_near] = i;
	matching[best_middle] = best_distant;
	matching[best_distant] = best_middle;
	(*red_n)--;
      }else if (gain_2>0){
	matching[i] = best_2;
        stack_free[free_p++] = matching[best_2];
        matching[matching[best_2]] = matching[best_2];
        matching[best_2] = i;
      } 
    } 
  }
  return 0;
}


/**
 * \brief MWM implementation
 *
 * MWM algorithm implementation, this is only an interface to the
 * N-cubed weighted matching Implementation of H. Gabow's Ph.D. thesis,
 * Stanford Univ. 1973 written by \author Edward Rothberg  7/85
 * \param coactrl coarsening control structure
 * \return coarsened graph
 */

GraphType* MWM(CoaCtrl* coactrl){

  GraphType *graph   = coactrl->g;
  int limit          = coactrl->npro;
  int randguide      = coactrl->randguide;
  int *guide         = coactrl->guide;  
  PartInfoType *info = coactrl->info;

  int i,j,nb_vtxs=0,ewgt=0,minvwgt=0;
  int jstart,jend,w;

  int  nvtxs  = graph->nvtxs;
  int* adjidx = graph->adjidx;
  int* adjncy = graph->adjncy;
  int* adjwgt = graph->adjwgt;
  int* vwgt   = graph->vwgt;
  int* map    = graph->map;
  int* match  = graph->match;

  Graph MWMgraph;
  int *Mate;
  int v;

  /************************************/
  /* init                             */
  /************************************/
  if(randguide == 1)
    randnum(nvtxs, guide);
  
  for (i=0; i<nvtxs; i++)
    map[i] = match[i] = -1;
  
  /************************************/
  /* translate GraphType to MWMgraph  */
  /************************************/
  MWMgraph = NewGraph(nvtxs);

  /* not to have to small weight using UPGMA */
  if(info->WType == IO_WEIGHT_UPGMA){
    for (i=0; i<nvtxs; i++)
      if(vwgt[i]<minvwgt || minvwgt == 0)
	  minvwgt = vwgt[i];
  }
      
  for (i=0; i<graph->nvtxs; i++) {
    v = guide[i];
    NLabel(MWMgraph,v+1) = v+1;
    Xcoord(MWMgraph,v+1) = 0;
    Ycoord(MWMgraph,v+1) = 0;

    jstart = adjidx[v];
    jend   = adjidx[v+1];

    for (j=jstart; j<jend; j++){
      w = adjncy[j];
      
      if(v<w){
	switch(info->WType){
	case IO_WEIGHT_UPGMA:
	  ewgt = 1000000 * adjwgt[j] / ((vwgt[w]-minvwgt + 1)*(vwgt[v] - minvwgt + 1));
	  break;
	case IO_WEIGHT_NORM:
	default:
	  ewgt = adjwgt[j];
	  break;
	}
	
	AddEdge(MWMgraph,v+1,w+1,ewgt);	
      }
    }
  }
  /************************************/
  /* compute matching                 */
  /************************************/
  Mate = Weighted_Match(limit,MWMgraph,1);

  /************************************/
  /* build graph->map graph->match    */
  /************************************/
  for (i=0; i<nvtxs; i++){
    v = guide[i];
    if(Mate[v+1]==0){
      map[v]   = nb_vtxs++;
      match[v] = v;
    }else if(Mate[v+1]<v+1){
      map[v]             = map[Mate[v+1]-1]=nb_vtxs++;
      match[v]           = Mate[v+1]-1;
      match[Mate[v+1]-1] = v;
    }    
  }  
  
  /************************************/
  /* contract the graph                */
  /************************************/
  graph = Contract(graph, guide, nb_vtxs);

  /*end*/
  return graph;
}
/**
 * \brief Modified Heavy Edge Matching implementation from \authir Wesley
 *
 * \param coactrl coarsening control structure
 * \return coarsened graph
 */
GraphType* MHem(CoaCtrl* coactrl){
  GraphType *graph   = coactrl->g;
  int limit          = coactrl->npro;
  int randguide      = coactrl->randguide;
  int *guide         = coactrl->guide;
  int *label         = coactrl->label;
  PartInfoType *info = coactrl->info;
  int maxvwgt        = coactrl->maxvwgt;

  int nvtxs          = graph->nvtxs;
  int* adjidx        = graph->adjidx;
  int* adjncy        = graph->adjncy;
  int* adjwgt        = graph->adjwgt;
  int* vwgt          = graph->vwgt;
  int* map           = graph->map;
  int* match         = graph->match;

  int i, j, jstart,jend, z, v, w, wstart, wend;
  double ewgt=0.0;
  int vleft; 		/*sum of the vertices which have not 
			  been matched*/
  int countwgt;		/*auxiliaryvariable holding the sum of edge 
			  weights to all adjacent vertices of a vertex*/
  int count;		/*auxiliary variables for counting number of 
			  vertices having the same maximum edge weight 
			  to a vertex*/
  double maxwgt = 0.0;	/*auxiliary variables recording the maximum 
			  edge weight to a vertex*/
  int nb_vtxs = 0;	/*number of vertices in the coarser graph */
  int s_label;		/*selected vertex */
  int s_wgt;		/*selected vertex sum of weight */
  
  /************************************/
  /* init                             */
  /************************************/
  if(randguide == 1)
    randnum(nvtxs, guide);
  
  for (i=0; i<nvtxs; i++)
    map[i] = match[i] = -1;
  vleft = nvtxs;

  /************************************/
  /* compute matching                 */
  /************************************/
  for (i=0; i<nvtxs; i++) {

    v = guide[i];
    	
    /*If the vertex v is still unmatched*/
    if (match[v] == -1){ 
	    
      maxwgt = 0.0;
      count = 0;

      jstart = adjidx[v];
      jend   = adjidx[v+1];
      
      /* add in label all adjacent vertices of v that have max edge weight*/
      for (j=jstart; j<jend; j++){
	w = adjncy[j];

	if (match[w] == -1 && vwgt[v] + vwgt[w] <= maxvwgt){ /*vertex not matched*/
	  switch(info->WType){
	  case IO_WEIGHT_UPGMA:
	    ewgt = (double)adjwgt[j] / (double)(vwgt[w] * vwgt[v]);
	    break;
	  case IO_WEIGHT_NORM:
	  default:
	    ewgt = graph->adjwgt[j];
	  }
	  
	  /* find edges of maximum weight */
	  if (maxwgt < ewgt){
	    count = 0;
	    maxwgt = ewgt;
	    label[count++] = w;
	  }else if (maxwgt == ewgt){
	    label[count++] = w;
	  }
	}
      }	

      /* MHEM heuristic : 
	 pick the adjacent vertex in label having maximum sum of weights
	 to all v's adjacent vertices. 
	 If there are more than one vertices meet this requirement, 
	 pick the first one */
      if(count == 0){		/* no possible match */
	s_label = v;
      }else if(count == 1){	/* only one choice */
	s_label = label[0];
      }else{	/* multiple choice */
	s_label = 0;
	s_wgt = -1;

	for (z=0; z<count; z++){
	  countwgt = 0 ;
		
	  jstart = adjidx[v];
	  jend   = adjidx[v+1];
	  for (j=jstart; j<jend; j++){
	    wstart = adjidx[label[z]];
	    wend   = adjidx[label[z] + 1];
	  
	    for (w=wstart; w<wend; w++){
	      if (adjncy[w] == adjncy[j]){
		countwgt += adjwgt[w];
		break;	/*break*/
	      }
	    }
	  }
	  if (countwgt > s_wgt){
	    s_label = label[z];
	    s_wgt = countwgt;
	  }
	}
      }
	
      /*Assign the values of map[][] and match[][]*/
      /*give a continues integer to be label of the new graph*/
      map[v]        = map[s_label] = nb_vtxs++;		
      match[v]      = s_label; 	
      match[s_label]= v; 	
      vleft -= 2;
      
      /*if the number of vertices is same as the number of processors, 
	match remaining vertices to themselves */
      if (nb_vtxs + vleft == limit){
	while (++i< nvtxs){
	  v = guide[i];
	  if (match[v] == -1){ 
	    map[v]   = nb_vtxs++;
	    match[v] = guide[i];
	  }
	}
      }
    }	
  }  

  /************************************/
  /* contract the graph               */
  /************************************/
  graph = Contract(graph, guide, nb_vtxs);

  /*end*/
  return graph;
}
/**
 * \brief Heavy Edge Matching implementation
 *
 * \param coactrl coarsening control structure
 * \return coarsened graph
 */
GraphType* Hem(CoaCtrl* coactrl){
  GraphType *graph   = coactrl->g;
  int limit          = coactrl->npro;
  int randguide      = coactrl->randguide;
  int *guide         = coactrl->guide;
  PartInfoType *info = coactrl->info;
  int maxvwgt        = coactrl->maxvwgt;

  int nvtxs          = graph->nvtxs;
  int* adjidx        = graph->adjidx;
  int* adjncy        = graph->adjncy;
  int* adjwgt        = graph->adjwgt;
  int* vwgt          = graph->vwgt;
  int* map           = graph->map;
  int* match         = graph->match;

  int i, j, jstart, jend, v, w;
  double ewgt=0.0;
  int vleft; 		/*sum of the vertices which have not 
			  been matched*/
  double maxwgt = 0.0;	/*auxiliary variables recording the maximum 
			  edge weight to a vertex*/
  int nb_vtxs = 0;	/*number of vertices in the coarser graph */
  int s_label;		/*selected vertex */
  
  /************************************/
  /* init                             */
  /************************************/
  if(randguide == 1)
    randnum(nvtxs, guide);
  
  for (i=0; i<nvtxs; i++)
    map[i] = match[i] = -1;
  vleft = nvtxs;

  /************************************/
  /* compute matching                 */
  /* HEM heuristic :                  */ 
  /* pick the adjacent vertex         */
  /* having maximum edge weight       */
  /************************************/
  for (i=0; i<nvtxs; i++) {

    v = guide[i];
	
    /*If the vertex v is still unmatched*/
    if (match[v] == -1){ 
	    
      maxwgt = 0.0;
      s_label = v;

      jstart = adjidx[v];
      jend   = adjidx[v+1];

      for (j=jstart; j<jend; j++){    
	w = adjncy[j];
	
	if (match[w] == -1 && vwgt[v] + vwgt[w] <= maxvwgt){ /*vertex not matched*/
	  switch(info->WType){
	  case IO_WEIGHT_UPGMA:
	    ewgt = (double)adjwgt[j] / (double)(vwgt[w] * vwgt[v]);
	    break;
	  case IO_WEIGHT_NORM:
	  default:
	    ewgt = adjwgt[j];
	  }
	  
	  /* find edge of maximum weight */
	  if (maxwgt < ewgt){
	    s_label = w;
	    maxwgt = ewgt;
	  }
	}
      }	

      /*Assign the values of map[][] and match[][]*/
      /*give a continues integer to be label of the new graph*/
      map[v]         = map[s_label] = nb_vtxs++;		
      match[v]       = s_label; 	
      match[s_label] = v; 	
      vleft -= 2;
      
      /*if the number of vertices is same as the number of processors, 
	match remaining vertices to themselves */
      if (nb_vtxs + vleft == limit){
	while (++i< nvtxs){
	  v = guide[i];
	  if (match[v] == -1){ 
	    map[v]   = nb_vtxs++;
	    match[v] = guide[i];
	  }
	}
      }
    }	
  }

  /************************************/
  /* contract the graph                */
  /************************************/
  graph = Contract(graph, guide, nb_vtxs);

  /*end*/
  return graph;
}
/**
 * \brief Random Matching implementation from \authir Wesley
 *
 * \warning this is not trully random, a trully random implementation 
 * would be too slow (see source)
 *
 * \param coactrl coarsening control structure
 * \return coarsened graph
 */
GraphType* Rm(CoaCtrl* coactrl){
  GraphType *graph   = coactrl->g;
  int limit          = coactrl->npro;
  int *guide         = coactrl->guide;
  /*int* guide2        = coactrl->label;*/
  int maxvwgt        = coactrl->maxvwgt;
  int randguide      = coactrl->randguide;

  int nvtxs          = graph->nvtxs;
  int* adjidx        = graph->adjidx;
  int* adjncy        = graph->adjncy;
  int* vwgt          = graph->vwgt;
  int* map           = graph->map;
  int* match         = graph->match;

  int i, j, jstart, jend, v, w;
  int vleft; 		/*sum of the vertices which have not 
			  been matched*/
  double maxwgt = 0.0;	/*auxiliary variables recording the maximum 
			  edge weight to a vertex*/
  int nb_vtxs = 0;	/*number of vertices in the coarser graph */
  int s_label;		/*selected vertex */
  
  /************************************/
  /* init                             */
  /************************************/
  if(randguide == 1)
    randnum(nvtxs, guide);
  
  for (i=0; i<nvtxs; i++)
    map[i] = match[i] = -1;
  vleft = nvtxs;

  /************************************/
  /* compute matching                 */
  /* RM heuristic :                   */ 
  /* pick a random adjacent vertex    */
  /************************************/
  for (i=0; i<nvtxs; i++) {

    v = guide[i];
    	
    /*If the vertex v is still unmatched*/
    if (match[v] == -1){ 
	    
      maxwgt = 0.0;
      s_label = v;
  
      jstart = adjidx[v];
      jend   = adjidx[v+1];

      /* guide for adjacent vertices */
      /*randnum(jend - jstart, guide2);*//* for trully random matching */
      for (j=jstart; j<jend; j++){     
	/*w = adjncy[jstart + guide2[j - jstart]];*//* for trully random matching */
	w = adjncy[j];

	if (match[w] == -1 && vwgt[v] + vwgt[w] <= maxvwgt){ /*vertex not matched*/
	  s_label = w;
	  break;
	}
      }
    
      /*Assign the values of map[][] and match[][]*/
      /*give a continues integer to be label of the new graph*/      

      map[v]         = map[s_label] = nb_vtxs++;		
      match[v]       = s_label; 	
      match[s_label] = v; 	
      vleft -= 2;
      
      /*if the number of vertices is same as the number of processors, 
	match remaining vertices to themselves */
      if (nb_vtxs + vleft == limit){
	while (++i< nvtxs){
	  v = guide[i];
	  if (match[v] == -1){ 
	    map[v]   = nb_vtxs++;
	    match[v] = guide[i];
	  }
	}
      }
    }
  }	
  

  /************************************/
  /* contract the graph                */
  /************************************/
  graph = Contract(graph, guide, nb_vtxs);

  return graph;
}

